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We revisit the paradigm of an ideal gas under isothermal conditions. A moving piston performs 
work on an ideal gas in a container that is strongly coupled to a heat reservoir. The thermal coupling 
is modelled by stochastic scattering at the boundaries. In contrast to recent studies of an adiabatic 
ideal gas with a piston [R.C. Lua and A.Y. Grosberg, J. Phys. Chem. B 109, 6805 (2005); I. Bena 
et al., Europhys. Lett. 71, 879 (2005)], container and piston stay in contact with the heat bath 
during the work process. Under this condition the heat reservoir as well as the system depend on 
the work parameter A and microscopic reversibility is broken for a moving piston. Our model is thus 
not included in the class of systems for which the non-equilibrium work theorem has been derived 
rigorously either by Hamiltonian [C. Jarzynski, J. Stat. Mech. P09005 (2004)] or stochastic methods 
[G.E. Crooks, J. Stat. Phys. 90, 1481 (1998)]. Nevertheless the validity of the non-equilibrium work 
theorem is confirmed both numerically for a wide range of parameter values and analytically in the 
limit of a very fast moving piston, i.e. in the far non-equilibrium regime. 
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I. INTRODUCTION 

Despite its successful application in numerous exper- 
imental and numerical studies (see e.g. the review ar- 
ticles [HE]), the validity of the non-equilibrium work 
theorem or 'Jarzynski relation' [5] still remains under 
discussion. Indeed this ongoing critique is mainly due 
to the surprising nature of the theorem: it states an ex- 
act equality that holds in situations arbitrarily far from 
equilibrium, under very general assumptions. More pre- 
cisely, it states that the free energy difference between 
two equilibrium states can be extracted from work mea- 
surements along irreversible trajectories connecting these 
two states. Therefore one can, in principle, obtain equi- 
librium information from a non-equilibrium experiment, 
which is of particular interest in chemical and biophysical 
applications. For example, the non-equilibrium work the- 
orem has been successfully applied to the stretching of a 
single protein |Jj . In this experiment the work performed 
by a single RNA molecule tethered between a solid sub- 
strate and a controllable cantilever in an aqueous salt 
solution is measured for slow (reversible) and fast (irre- 
versible) stretching. The free energy difference between 
its folded and unfolded conformations is obtained from 
the reversible process using ordinary equilibrium ther- 
modynamics. On the other hand applying Jarzynski's 
relation to the work values obtained from the irreversible 
process also reproduces this result within experimental 
errors, thus confirming the theorem. 

However it has been questioned, in Ref. [5], whether 
this experiment indeed creates a non-equilibrium situa- 
tion. It is argued that a slow or fast work process does 
not necessarily guarantee its reversibility or irreversibil- 
ity. Rather the work rate has to be compared with the 
strength of the coupling (rate of heat transfer) between 
the system and its thermal environment. If the work 
rate is apparently large, but still smaller than the rate 



of heat transfer, the system is essentially maintained in 
an equilibrium state. This is claimed [5] to be the case 
in the protein stretching experiment, since the surround- 
ing liquid allows for rapid thcrmalization. Under such 
conditions the theorem is expected to hold trivially. 

The above discussion highlights the importance of 
properly assessing the thermostating process between 
the system and the heat reservoir. The purpose of the 
present paper is to investigate the Jarzynski relation for 
the most simple thermodynamic system under isother- 
mal and non-equilibrium conditions. In a gedankenex- 
pcrimcnt an ideal gas is isothcrmally expanded in a heat- 
conducting container by pulling a piston at different ve- 
locities. Work is performed when the gas particles hit 
the piston during its movement. Similar ideal gas mod- 
els have been investigated in [5J |7] , but under adiabatic 
conditions, i.e. without considering heat transfer during 
the work process. The extension to an isothermal situ- 
ation provides important further insight into Jarzynski's 
relation. 

The remainder of this paper is organized as follows. In 
the next section we present a brief outline of the non- 
equilibrium work theorem for a system with strong ther- 
mal coupling. In Sec. |III| the isothermal ideal gas model 
is introduced, which allows for an analytical formulation 
of the non-equilibrium work theorem. We present the re- 
sults of a numerical study of this model and revisit the 
adiabatic expansion of the ideal gas in Sec. |IV| Finally 
we conclude with a summary of the main points and a 
brief outlook. 



II. THE NON-EQUILIBRIUM WORK 
THEOREM 

The non-equilibrium work theorem can be formulated 
in the following way. The system of interest is prepared 
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in an initial state of equilibrium while in contact with a 
heat reservoir at temperature T. By changing an exter- 
nal parameter A — the work parameter — according to 
a fixed protocol from A to B the system is subjected to 
a thermodynamic process at the end of which it reaches 
a final state not necessarily in equilibrium. This pro- 
cess can possibly drive the system arbitrarily far away 
from equilibrium while performing a certain amount of 
work, W. During the work process the system may or 
may not stay in contact with the thermal environment. 
Upon reaching the final parameter value A = B the sys- 
tem relaxes to equilibrium by exchanging heat with the 
reservoir but it is assumed that no further work is per- 
formed. If we repeat this process following the same pro- 
tocol infinitely many times we obtain a distribution of 
work values p(W) due to the stochastic nature of the ini- 
tial equilibrium state, which is sampled from a canonical 
distribution. The Jarzynski relation then states a strong 
constraint on this work distribution |3J: 

(e-W) = JdW p(W) e-? w = e-? AF . (1) 

The average over the exponentiated work values equals 
the exponential of the free energy difference between the 
initial and final equilibrium states of the system, whether 
or not the final equilibrium state is actually realized in 
the disturbed system. j3 is the inverse temperature 1/(3 — 
ksT and AF is the ratio of the equilibrium partition 
functions: 



A. Hamiltonian derivation 

In Jarzynski's original derivation [3] almost a decade 
ago, the specific assumption was made that the coupling 
between system and heat reservoir is sufficiently small to 
neglect the interaction term in the Hamiltonian. Under 
this condition the work process is effectively assumed to 
be adiabatic. Only recently a derivation has been pre- 
sented that does not rely on the weak coupling assump- 
tion [5]. The starting point is the Hamiltonian 

H(T; A) = H(x; A) + H E (y) + h int (x, y), (3) 

for the system and thermal environment. H (x; A) denotes 
the Hamiltonian of the system, H E (y) is the Hamilto- 
nian of the heat reservoir and hi n t(x,y) is the interac- 
tion Hamiltonian. Here, x refers to a point in the phase 
space of the system only, likewise y for the heat reser- 
voir. We denote a point in the combined phase space by 
r = (x,y). It is important to note that the dependence 
on the work parameter enters only via the Hamiltonian 
of the system, while the heat reservoir is assumed to be 
A-indcpcndcnt. The combined system-plus-reservoir is 
now subject to an adiabatic work process A(i), such that 
the non-equilibrium work theorem can be applied in its 
original form. For finite h int (x, y) the equilibrium distri- 
bution of the system has to be described by the modi- 
fied Boltzmann-factor ps{x; A) oc exp[— (3H*(x; A)] where 
H*(x; A) is called a potential of mean force [21 E]: 



AF = F B -F A = -p- l \n^. (2) 

& A 

It should be noted that Fb corresponds to the actual 
free energy of the final state only if the work parameter 
is finally held fixed at A = B until the system has ther- 
malized with the reservoir. Without thermalization the 
system may well be in a state out of equilibrium such 
that no free energy can be defined. 

The relation, Eq. ([lj, holds irrespective of the partic- 
ular character of the work process and is valid beyond 
the linear response regime. For a fast switching of A one 
may perform a non-equilibrium process and still obtain 
the equilibrium free energy difference by evaluating the 
exponential work average. Thus the Jarzynski relation 
is one of the few exact results applicable far from equi- 
librium. Since (exp(x)) > exp((x)), Eq. implies the 
second law of thermodynamics formulated for work and 
free energy, (W) > AF. The equality (W) = AF is only 
true for a reversible (quasistatic) process. 

The Jarzynski relation has been derived both for deter- 
ministic Hamiltonian dynamics [3, 8 and for stochastic 
dynamics [9] . Without resorting to the full proof we shall 
make a few comments on the derivation of Eq. for a 
system with strong thermal coupling. 



H*(x; A) 

tt( o-i, fdyexp[-/3(H E (y) + h int (x,y))} 

= H(x; A - (3 In J - j- Wa~T\\ 

J &yexv[-f3H E {y)\ 

(4) 

With this consideration, the left hand side of Eq. 
is reduced to the ratio of the partition functions of the 
system only, resulting in the same result as in the case of 
weak coupling. 

One should note a subtlety that applies to systems with 
rigid boundaries whose position varies with the control 
parameter X(t) (which is the case e.g. for the free ex- 
pansion of an ideal gas in a box, or expansion against a 
piston). Such boundaries have to be regarded as poten- 
tials in the Hamiltonian of the system. If one imposes 
instead time-dependent constraints, the resulting Hamil- 
tonian evolution does not conserve phase space volume 
thus leading to an apparent violation of Jarzynski's rela- 
tion. The correct procedure assumes a potential strength 
depending on a parameter e which in the limit e — > oo 
becomes a rigid boundary. As a consequence evaluating 
the exponential work average in Eq. (fTl requires correct 
ordering of the thermodynamic limit (number of repeti- 
tions) before the e limit [T5] . 
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B. Derivation for stochastic microscopic reversible 
dynamics 

In [3] the non-equilibrium work theorem was derived 
under the assumption of stochastic Markovian microscop- 
ically reversible dynamics. The crucial condition of mi- 
croscopic reversibility is formulated as 9 : 



P(x(t)\X(t)) 

p(x(-t)\M-t)) 



= exp[-/?Q(z(t),A(t))]. 



(5) 



Here, x(t) is a particular trajectory of the system when it 
is subject to the external work process A(t). P(x(t)\X(t)) 
is the probability of following this path during the work 
process and P(x(— t)|A(— t)) is the probability of the cor- 
responding time-reversed path during the time-reverse of 
the process. Q denotes the heat transferred from the 
heat reservoir to the system during this process; it is a 
functional of the path, with the property Q(x(t), A(i)) = 
Q(x(—t),X(—t)) under time reversal. Equation d5J) is 
claimed to hold arbitrarily far from equilibrium and al- 
lows for a simple proof of Jarzynski's relation as discussed 
in [5J [TH] . In the following we discuss the derivation of 
Eq. (§ more thoroughly. 

Microscopic reversibility is usually expressed by the 
condition of detailed balance for the transitions between 
states a and b [TT] : 



uj{b — ► a) 



(6) 



where the 'states' a, b refer to discrete volume elements 
in the phase space of the system and Lo(a — ► 6), u>(b — ► a) 
are the corresponding transition rates. It is argued in [9] 
that Eq. ([5]) follows directly by discretizing the path into 
single time steps, which each obey detailed balance. More 
precisely, the evolution of the system x(t) = {xq, ■■■,Xt} 
is considered for a fixed sequence of the work param- 
eter X(t) = {Ai,...,A t } such that the single time steps 
can be decoupled into two substeps. First the control 
parameter is changed, Xi — ► Aj+i, performing a certain 
amount of work, and then the system evolves for fixed 
Aj+i, Xi — > Xi + \, exchanging heat with the thermal envi- 
ronment. As a consequence of Markovian dynamics the 
probability P(x(t)\X) of following a path through phase 
space under the work process can be expressed as a prod- 
uct of transition probabilities for the discretized quanti- 
ties {x , ...,x t } and {Ai, A t } 0: 



t-i 



P(x(t)\X(t)) = Y[P{xi -> ar i+ i|A i+1 ). 



(7) 



Then the ratio of probabilities of a forward path and its 
corresponding time-reversed path is 



P(x(-t)\X(-t)) 



n 



P(xj -> x l+1 \X i+1 ) 
P(xi <- ac i+1 |Ai+i) 



-0Q(x(t),\(t)) 



(8) 



In the second line it is assumed that, for each time step, 
a 'detailed balance like' condition holds for the ratios of 
forward and reverse probabilities analogous to Eq. @ for 
fixed values of the work parameter X. The third line is 
due to the first law of thermodynamics: the difference in 
energy between two successive states is completely sup- 
plied by the heat bath if no work can be performed (i.e. 
for constant A), AQ = E(xi+i, A,+i) — E[xi, Aj+i). Sim- 
ilarly the work performed by the system originates only 
from a change in A: AW — E(xi, Aj+i) — E(xi, Xi). This 
decoupling is crucial for the derivation of Eq. (J5|. By 
considering only transitions x% — > Xi+\ for constant work 
parameter {Ai, A t } in Eq. (|8]), the evolution of the sys- 
tem is effectively reduced to a sequence of static states 
that obey detailed balance. Crucial to this is the assump- 
tion that the transition rates are independent of the rate 
A at which the system is disturbed. One might there- 
fore question whether Eq. ([5]) is indeed valid away from 
equilibrium, where detailed balance is not generally ex- 
pected to hold (T2j [13]. If detailed balance is violated 
under particular conditions, then Eq. (|5]) also fails. It 
is thus reasonable to state that the derivation of Eq. ^ 
is correct under the given assumptions, but that these 
assumptions do not properly take into account a non- 
equilibrium evolution of the system. 

In the next section we will further discuss micro- 
scopic reversibility with regard to the isothermal ideal 
gas model. 



III. THE IDEAL GAS WITH A PISTON: AN 
ISOTHERMAL MODEL 



We consider a one-dimensional classical non- 
interacting ideal gas in a container with a moving 
piston (as shown in Fig. [I] which resembles the system 
treated in reference [5]). Both the end wall of the con- 
tainer and the piston are connected to a heat reservoir 
which keeps the gas at constant temperature. This heat 
reservoir is assumed to be in thermal equilibrium such 
that its degrees of freedom are distributed according to a 
canonical Boltzmann distribution. Interactions between 
the system and its thermal environment are modeled by 
stochastic scattering at the boundaries: when the gas 
particle reaches either side of the container a completely 
inelastic collision takes place and the particle loses all 
its kinetic energy. It receives a new stochastic velocity 
which is sampled from the probability distribution of 
the heat reservoir, independent of its former velocity. 
We refer to this situation as a strong thermal coupling 
between system and environment. 

Due to flux conservation, the probability distribution 
of the new particle velocity v out after a particular col- 
lision with the fixed (left-hand) boundary (proportional 



FIG. 1: The ideal gas confined in a container with a piston. 
The initial position and velocity of a single gas particle are 
denoted x and «i, and the piston has velocity v p . The (one- 
dimensional) gas is initially confined to a length L. 



to the flux of particles leaving the boundary) therefore 
takes the form [16] : 

M^.) = !gfU(-^). (9) 

Here Vi always denotes the modulus Vi = Under 
equilibrium conditions, this boundary condition yields 
the Boltzmann distribution for the velocities of particles 
in the container volume. A similar expression applies 
to the distribution of velocities assigned at the moving 
boundary (the piston), but with the important distinc- 
tion that this is the distribution of out-going velocities in 
the frame of the piston. Movement of the piston there- 
fore results in a non-zero mean streaming velocity in the 
laboratory frame. 

In contrast to the commonly used Gaussian or Nose- 
Hoover thermostating schemes, this stochastic boundary 
thermostat is non-deterministic and non-time-reversible 
(in the lab frame) . Nevertheless it provides a valid phys- 
ical model of the heat bath interaction. Furthermore, 
since no potential acts on the ideal gas particles, their 
energy is purely kinetic and completely determined by 
the canonical probability distribution of the heat reser- 
voir. 

An important property of the heat reservoir in this 
model is its dependence on the work parameter A, which 
is more precisely a A-dependence of its center of mass. In 
the counter-intuitive context of the Jarzynski equation 
this should not be considered trivial. As has been men- 
tioned in the previous section, the derivation of the non- 
equilibrium work theorem assumes A-dependence only for 
the Hamiltonian of the system, not for the heat reser- 
voir. Here, that assumption is violated by the moving 
piston, which, by definition, changes its position as a 
function of A and, at the same time, thermostats the 
system. Although a rigorous and general treatment of 
this issue would require a Hamiltonian description of the 
heat reservoir, the stochastic model provides important 
insight. 

The isothermal ideal gas is subjected to the following 



FIG. 2: Position-time-diagram for the work process. The 
particle performs work by hitting the piston during the time of 
its movement r. Example trajectory for a positive (negative) 
initial velocity is denoted by a solid (dotted) line. Note the 
difference from the corresponding diagram in [BJ. Here the 
particle trajectories have varying slopes due to the isothermal 
boundaries. 



thermodynamic process (see Fig. [2]). In the initial state 
the gas is in equilibrium with the heat bath, confined to 
the (one-dimensional) 'volume' L by a fixed position of 
the piston. An individual gas particle samples its veloc- 
ity from a Maxwell-Boltzmann distribution pmb{v) and 
its position from the uniform distribution l/L. Then the 
external work process begins by moving the piston out- 
wards at constant speed v p for a time period r. Work 
is performed when gas particles collide with the piston 
during its movement. (The retreating piston then does a 
negative amount of work on the gas.) When the piston 
stops the gas thermalizes to the final equilibrium state at 
volume L + v p T. 

With regard to the previous discussion we identify the 
work parameter A as the position of the piston and the 
work rate as dependent on A = v p . Since the speed of the 
heat-transfer mechanism remains fixed in this model, the 
piston velocity determines the reversibility or irreversibil- 
ity of the process. In the limit v p — > with r — > oo, we 
perform a quasistatic reversible expansion of the gas. In 
the converse case of large piston speed and short t, the 
bulk of the gas remains in the initial part of the con- 
tainer after the piston has stopped, although the volume 
is extended. Subsequent equilibration is only completed 
sometime after the full volume has been explored by the 
gas. 



A. Microscopic reversibility 

It has been pointed out in Sec. |IIB| that the validity 
of the non-equilibrium work theorem for stochastic dy- 
namics in [5J [TO] relies on the condition of microscopic re- 
versibility. The isothermal ideal gas model offers the op- 
portunity to check Crooks' assumption (Eq. Q) against 
a physically motivated, realistic case. In order to deter- 
mine the transition rates of Eq. (|6J) we have to discretize 



5 



the phase space of the system into intervals of size Ax, 
Av. As the internal energy is only changed by inelas- 
tic collisions with the boundaries, the discussion can be 
reduced to a single collision event in the vicinity of the 
piston. 

First we consider the equilibrium case, where the work 
parameter (i.e. the position of the piston) is fixed at L. 
The particle is in state a if it occupies the phase space 
element [v a , v a + Av] x [L — Ax, L] and in state b if it 
is found in [v b , Vb + Av] x [L — Ax, L]. We make Av in- 
finitesimally small. Let a and b refer to the same position 
interval (but different velocities) before and after the col- 
lision. Obviously, in order to make the transition a — ► b 
by a bounce against the piston, v a and v& have opposite 
directions. 

A transition rate u>(a — ► b) = p(b, At\a,0) / At is given 
in terms of a conditional probability p(b, At\a, 0) for the 
system to be found in state b at time At, given that it 
was in state a at initial time to = 0. In turn this condi- 
tional probability is determined by the joint probability 
p(b, At; a, 0) of finding the particle in state a at time 
and in state b some time At later: 



p{b, At\a,0) 



p(b, At;a,0) 

pM) 



(10) 



The phase space probability distribution of an ideal gas 
at equilibrium is p e? (x, v) = 1/L pmb(v), thus the proba- 
bility p{a, 0) is simply given by p(a, 0) = p eq (x, v a )AxAv. 
On the other hand, the calculation of the joint proba- 
bility p(b, At; a, 0) involves the following considerations. 
Initially the particle occupies state a. Relative to the 
piston its position is x = x — L with modulus x <E [0, Ax] 
and its velocity is v a . After the bounce the new particle 
velocity is sampled from p B (Eq. ([9])). Given that v a and 
Vb are fixed, the transition a — > b only occurs if x fulfils 
two conditions. First, the collision time t c = x/v a < At, 
since the collision has to take place within At. Second, 
when the particle travels at V(, after the collision, it is 
found within [0, Ax] at time At only if Ax > (At — t c )vb- 
Consequently the joint probability p(b, At; a, 0) is deter- 
mined as: 

p(b,At;a,0) = (O(x)O(Ax - x)0(At - t' c ) 
x©(Ax- (At -t' c K) 
xS(v a -v' a )6(v b -v' b )). (11) 

The average is taken with respect to the stochastic vari- 
ables x, v' a , v' b , which are sampled from the distributions 
p eq (L - x,v' a ) and p B (v' b ): 

p(b,At;a,0) = p B (v b )Av dx ©(At - x/v a ) 

Jo 

xO(Ax - (At - x/v a )v b ) 

xp eq (L - x,v a )Av. (12) 

Since the particle position is uniformly distributed, 
p eq {L — x, v a ) can be taken out of the average and finally 



cancels when the conditional probability Eq. (10) is con- 
sidered. One is essentially left with an integral over a 
product of two theta functions, whose result is presented 
in appendix [B] We obtain for the conditional probability 
p(6,Af|o,0): 



p(b,At\a,0) = A{v a ,v b )e~ v "' 2 , 



(13) 



where A(v a ,v b ) is a function symmetric in v a , Vb (k B T 
has been set to unity). The calculation of the condi- 
tional probability p{a, At|6,0) for the transition b — > a 
under time reversal follows the same steps and yields 
Eq. (13 1 with v a , v b exchanged. It then follows that the 



ratio of transition rates between states a and b is equal 
to the Boltzmann factor of their energy difference. We 
are therefore reassured that, in the absence of external 
work, the model respects the detailed balance condition, 
Eq. {6|. 

In the non-equilibrium case, the movement of the pis- 
ton has to be taken into account. In the context of the 
discussion in Sec. |IIB[ the derivation of Crooks' micro- 
scopic reversibility condition, Eq. (|5j), is tantamount to 
fixing the work parameter (the piston position) at a suc- 
cession of different values L(t) = {L + v p t\, ...,L + v p r}, 
but disregarding the momentum exchange resulting from 
the movement. With that omission, the detailed balance 
condition follows trivially for each fixed piston position 
L+Vpti since we can simply repeat the calculations above 
with the new position L = L + v p ti . Equation ^ would 
then follow. 

However, the transition rates in the non-equilibrium 
case have to include the dynamic change of the work 
parameter. We will show that this leads to a violation of 
the balance condition used in the derivation of Eq. ([5j. 
Let us assume that the ideal gas is in a macroscopic non- 
equilibrium state due to the movement of the piston. The 
probability distribution of the particle is thus given by 
a function p{x,v), whose precise form is unknown. We 
consider the transition a — > b during time step At due 
to a collision event in the vicinity of the piston at initial 
time to- The position of the piston is L' = L + v p t and 
we can define state a similar to the static case above as 
the volume clement in front of the piston at to. Likewise 
state b is defined as the volume element in front of the 
piston at time to + At: 



a = [v a ,v a + Av] x [L 1 — Ax, L'], 
b = [v b , v b + Av] x [L 1 



v p At - Ax,L' - 



v p At](U) 



As before we are interested in the conditional probabil- 
ity p Vp (b,t + At\a,t ). For this purpose we transform 
to coordinates in the frame of the piston (all quantities 
referring to this frame will be denoted by an overbar): 



(15) 



As well as x' = x — (L + v p t). If we consider states a 
and b in piston coordinates: o = [v a , v a + Av] x [0, Ax], 
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b = [vb, Vb + Av] x [0, Aa;], we realize that they refer to the 
same situation as in the static case treated above. Conse- 
quently the joint probability p Vp (6, t + At; a, t ) assumes 
the form of Eq. [TTj The average is now taken with respect 
to x, v' a , v' b , which are sampled from p(JJ — x, v' a + v p ) 
and pB(v'b) (the latter distribution remains unchanged 
because v' b is the velocity obtained from the moving pis- 
ton). This average can be calculated as before, assuming 
that the non-equilibrium distribution p(JJ — x,v' a + v p ) 
becomes uniform with respect to x within the interval 
[0, Ax] for small Ax. As a consequence p(L' — x,v' a + v p ) 
can again be taken out of the average and cancels with 
the identical term in p(a,to). The conditional probabil- 
ity PvM>, t + At\a, t ) finally assumes the same form as 
Eq. (13 1. For the reverse transition b — > a under time 



reversal, the same considerations lead to the result as in 
the static case. We therefore find, that in the frame of the 
piston, the detailed balance condition is indeed fulfilled: 

Pv P (b,t + At\a,t ) _ A(v a ,v b ) _ c 2 /2+c 2 /2 



p Vp (a,t \b,t + At) A(v b ,v a ) 

_ -(E b -E a ) 



(16) 



However, the situation is different in the laboratory 
frame. Using the transformation Eqs. (15) in order to 



transform Eq. (16) to the lab frame, we sec that detailed 



balance is violated: 

PvJb,t + At\a,t ) 



Pv p (a,t \b,t + At) 



e -(v b + Vp y/2 
e -K-« P ) 2 /2' 



(17) 



Consequently Eq. ([5| does not hold either. As a result 
we observe that, in the laboratory frame, the isothermal 
ideal gas model with a moving piston violates Crooks' 
microscopic reversibility condition. 



B. The exponential work average 

In order to verify the non-equilibrium work theorem 
for the ideal gas with stochastic boundary conditions, the 
main task is to evaluate the exponential work average of 
the isothermal expansion process. The free energy dif- 
ference on the other hand can be calculated from simple 
thermodynamic considerations: 



AF = ln^ = iVlnf— — 
Zb \L + v p t 



(18) 



Throughout the calculations we set the temperature pa- 
rameter (3—1 without loss of generality. In the case 
of a non-interacting gas the partition functions Za (ini- 
tial state) and Zg (final equilibrated state) both factorise 
and one can effectively reduce the calculations to a sin- 
gle particle N = 1. This gas molecule performs work if it 
can hit the piston during its movement, i.e. during the 
time period r. Since it obtains a new randomly chosen 
velocity upon reaching either the wall or the piston, the 



time tk for k bounces to occur depends on all realized ve- 
locities vi,...,Vk such that tk = tk({v%, Vk},x; L, v p ). 
Throughout this article we shall use the convention that 
the variables Vi always refer to the velocity of the particle 
in the reference frame of the particular wall, where this 
velocity has been obtained. Also all calculations are per- 
formed with the modulus of the velocities. Since the new 
velocity of the particle after collision with the piston is 
thus given in the frame of the piston, the following recur- 
sion relation holds, as can easily be verified by inspecting 
the position-time diagram (Fig.[2|: 



Vk 



-tk-1 



(19) 



As a consequence we have to distinguish between a 
positive or negative initial velocity v± since, in each case, 
different bounces contribute to the work average. For a 
positive initial velocity the odd-numbered bounces yield 
the work contribution and the time for the first bounce 
is tf = (L — x)/{y\ — v p ); for a negative sign the even 
bounces contribute with t± = xjv\. Hence the exponen- 
tial work average can be written as: 

(e~ w ) = <e°) + + (e°>_ + {e~ w ) + + (e~ w )_ , (20) 

where +/— refers to a positive or negative initial veloc- 
ity and the average is to be taken over the uniformly 
distributed initial position x and all velocities vi,...,Vk- 
The zero-work contributions to the exponential work av- 
erage are then determined according to: 



1 f f°° 
(e°) , = - / dx / dv 1 pmb(vi) 
L Jo Jo 



[e(t+-r)e(*+) + e(-i+)] 



i 



(e°)_ = - dx dv 1 pmb(v\) / dv 2 Pb{v 2 ) 
L Jo Jo Jo 

(21) 



x [Q(q ~r)Q(q) + Q(-q)] 



Here <d(x) denotes the Heaviside theta function. Recall 
that a negative amount of work is performed on the gas 
by the piston, since we consider an expanding volume. 
The total work W is given by summation of the contri- 
butions due to the individual bounces against the piston. 
The single bounce contributions, 



= -Ap • v p = -{vi + v i+ i)v p + vf } 



(22) 



are determined by the momentum transfer Ap = Vj — 
(vj+i + v p ), where we have set the particle mass to unity 
without loss of generality. Here, Vj+i is always the new 
velocity of the particle after collision with the piston and 
is therefore given in the frame of the piston. The shift 
— v p takes this into account when considering the mo- 
mentum transfer in the laboratory fram e, le ading to the 



term +vi in the work contribution Eq. ( 22 1 . The Wi are 



statistically independent random variables 



7 



The average of the exponential work can be expressed 
as a series in the number of bounces n with the piston: 



1 °° pL p°° /*°° 

(e- W ) + = tE/ ^ **•"/ d 
n=l Jo J ° J ° 

P+ 2n - 1; r) e" t " 2i " 1 

(23) 

We introduce the joint probability distribution func- 
tion P + ({vx, V2n+i}, 2n — 1; r) which determines the 
probability of a particular realization of velocities 
{«!, ...,U2n+i} and °f exactly n bounces with the piston 
resulting, within the time period r: 

P + ({vi,..., v 2 „+i}, 2n - l;r) 
= Pmb(vi)pb(v2).-Pb(v2ti+i)[&(t - tan-l) x 

e(t+)...e(<+ _j - e(T - *J l+ i)e(t+)...e(t+ +1 )] 

(24) 

An analogous expression holds for {e~ w )_: 

^ oo „2j ~oo poo 

(e~~ W )_ = yY] / da: / dw i--- / d "2»+2 x 

P_ ({«!, U2 „ +2 }, 2n; r) e" EIU — (25) 
with the joint probability distribution 
P_ ({«!,..., i^ n+2 },2n;r) 

= PMb{vi)pb(V2):. PB(V2n+2)[@(T - fcfj X 

e(<r)...e(t 2 - n ) - 9(t - t 2 ; +2 )9(fr)...e(f 2 ; +2 )] 

(26) 

The difficulties in evaluating the exponential work av- 
erage originate primarily from the integration over the 
theta function 0(r — ife({t>i, a;)) where is given 

by the recursion relation ( 19 1 . Therefore we resort to 



V2n+l x 



a numerical investigation of the average, Eq. (20 1, in 
Sec. |IV| and in Sec. |III C| below to an analytical but ap- 
proximate evaluation, which tends to the exact answer in 
a well controlled limit. 



C. The limit of a fast moving piston 

The assumption that at most one bounce takes place 
between particle and piston yields an approximation that 
is analytically tractable and becomes exact in the limit 
of a fast moving piston while the volume extension is 
small compared with the original volume. In this case, 
L 3> v p T 3> t where velocities are measured in units 
of the thermal velocity since = m = 1. In Ref. [B] the 
one-bounce approximation validated the non-equilibrium 
work theorem for the adiabatic ideal gas expansion by 
considering the n = 1 approximation of the work distri- 
bution. Here we establish this result by calculating the 



one-bounce approximation of the exponential work aver- 
age Eq. pi)] ). According to the Jarzynski relation this 
average should yield the exponential of the free energy 
difference Eq. (18 1: 



e- AF = 1 



VpT 

L ' 



(27) 



The n — 1 approximation of the exponential work aver- 
age consists of four contributions: 

(e-^> n=1 = <e°) + + <e°>_ + <e-*>. + <^ 2 >_ -(28) 



Here (e°) and (e°)_ are given as in (21 1. The two 
one-bounce contributions, from Eqs. (23 1 and (25), are: 



L 



dx 



dvi pmb{v\) I dv 2 p B (v 2 ) 



x0(r-t+)0(t+) 



In the appendix we derive the following results in the 
limit L 3> v p T 3> t: 



dx dvip M B(vi) / dv 2 p B (v 2 ) 
L Jo Jo Jo 

poo 

x/ dv 3PB (v 3 )e(T-q)e(q)e~ W2 (29) 



<e°> + + <e >- - 1, 
{e~ W2 )_ -> 0. 



(30) 



In this section we calculate the dominant contribution 
(e~ Wl ) + . This average takes the explicit form: 



. dv 2 v 2 e - {v2 - v " )2/2 / dx 
2ttL Jo Jo 

L-x 
v\ - v p 



V v 1 -v p ) 



(31) 

The integrations over the v\ and v 2 variables are inde- 
pendent. For ^ > 1 the integral over v 2 yields ^/2ttv p . 
Making the substitution v = v± — v p and changing the 
argument of the theta function, we can calculate the re- 
maining integral in a straightforward way: 



I! 



dv e~ v /2 / dx e(x - (L - vt)) 
o Jo 

V -f [ L/T dvve-* 2 /* 
L Jo 

VpT 



r„ I dve~ v ' /2 

L/r 



(32) 



The last limit holds for L ^> v p t such that, in com- 



bination with Eq. (|30j) , we obtain the expected result 
Eq. Q. Fig. [3] shows a plot of (e~ w ) n=l for L = 1, 
v p = 1 together with Eq. (|27|) . According to our nu- 
merical results presented below a piston velocity in this 
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FIG. 3: Comparison of the n = 1 approximation Eq. (281 
(dotted line) with the exponential free energy difference 
1 + v p t/L (solid line) for v p — 1. The triangles show the 
numerical results for the same parameter values. 



regime drives the system sufficiently out of equilibrium. 
One notes that here the one-bounce approximation repro- 
duces Jarzynski's relation for times up to t rs 1, so that 
this approximation actually holds in a broader regime 
than the limits considered above. For larger times r the 
single bounce is no longer sufficient to approximate the 
averaged exponential and higher order terms have to be 
taken into account. 



IV. NUMERICAL METHOD 

We have performed a numerical study of the isothermal 
ideal gas model in order to evaluate the exponential work 
average, Eq. (20), for arbitrary numbers of bounces n. 



This thermostated 'Molecular Dynamics' simulation es- 
sentially consists of picking random numbers from the ve- 
locity distributions pmb{v\), Ps(wj) and calculating the 
bouncing times 4* (Eq. (|l~9|)). When the velocity config- 
uration allows the particle to hit the piston one or more 
times within t, the corresponding work values (Eq. ( 22 1 ) 
are recorded. 

A general problem of the applicability of the non- 
equilibrium work theorem is the convergence of the aver- 
age. Since the exponential exp[— (3W\ emphasizes small 
work values, one effectively has to sample the far left 
tail of the work distribution. In our model, the work 
performed on the system is always negative and the av- 
erage is dominated by those events that lead to a large 
negative work contribution. A general discussion of the 
convergence problem can be found in [2] with particular 
focus on the ideal gas and piston, but for an adiabatic 
work process (clastic collisions). It was shown that, in 
this particular case, the number of realizations needed in 
order to sample the dominant part of the average grows 



FIG. 4: Reversible isothermal expansion of the ideal gas for 



exponentially with the system size [B] or, more gener- 
ally, this number is proportional to the exponential of 
the averaged work that is dissipated during the reverse 
process [TJ. It is not obvious whether these results can 
be directly applied to the isothermal ideal gas model un- 
der consideration. In contrast to the adiabatic case, the 
occurrence of many collisions in a particular realization 
does not necessarily imply a large work contribution (a 
dominant event), as the work depends not only on the in- 
coming velocity but also on the statistically independent 
outgoing velocity. 

The numerical results below have been obtained with 
typically 10 6 realizations per data point, which yields 
such excellent convergence of the exponential average 
that error bars have been omitted in the figures. All 
less relevant parameters are set to unity for simplicity 
and units are non-dimensional. This includes the length 
of the initial volume, L = 1, and the inverse temperature, 
(3 = 1, which sets the width of pu b and ps ■ Accordingly 
the thermal velocity is unity as well so, for piston speeds 
v p > 1 , we are in a regime where the tail of the initial ve- 
locity distribution contributes the work. The plots show 
the average exponential work for a given v p when t is 
varied from zero up to the extended volume v p t = L 
(i.e. the volume is doubled). 



A. Numerical results 

We report the following results. In the limit of a very 
slow (quasistatic) expansion we obtain the isothermal 
free energy difference, Eq. (18), from both the work av- 
erage (W), in accordance with the Second Law for a re- 



versible process, and the exponential work average (see 
Fig. [Ij). The quasistatic regime is found at v p < 0.001. 
If we pull the piston at a higher speed the work average 
deviates noticeably from the free energy difference, indi- 
cating the onset of dissipation. The dissipated work is 
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FIG. 5: Irreversible isothermal expansion of the ideal gas 
for v p = 1. One notices the dissipated work as the difference 
between the average work and AF. 



FIG. 6: Work distribution of the isothermal expansion pro- 



Wd — (W) — AF and we are effectively performing an 
irreversible non-equilibrium experiment. On the other 
hand the negative logarithm of the exponential work av- 
erage still agrees with the isothermal free energy differ- 
ence as predicted by the non-equilibrium work theorem 
(see Fig. [5]). This is the main result of our numerical 
investigation: that the Jarzynski relation holds in the 
non-equilibrium regime of our model (y p > 0.001) de- 
spite the fact that the model does not belong to the class 
of systems for which the theorem was derived, but has a 
more physically motivated coupling to its heat bath. The 
simulation shows excellent convergence up to v p ~ 1 . For 
higher piston speeds it becomes increasingly difficult to 
sample the tails of the velocity distributions. 

In Fig. [6j we present the full distribution of work val- 
ues determined numerically for v p = 0.1. The distri- 
bution exhibits a multi-peak structure, demonstrating 
that, in imposing only one constraint on the distribu- 
tion, the non-equilibrium work theorem does not confine 



it to adopt a simple shape. 



B. The adiabatic piston model revisited 

An ideal gas with a piston was previously investigated 
by Lua and Grosberg |6j for the case of adiabatic ex- 
pansion, i.e. perfectly elastic collisions at the boundaries. 
We can reproduce their adiabatic model by considering 
elastic, energy-conserving (and therefore deterministic) 
collisions instead of completely inelastic, stochastic ones. 
Consequently the probability distribution, Eq. is sub- 
stituted by 



PB(v ut) = S(v out - fan ~ v p )), 



(33) 



when the incoming velocity is v in . The shift —v p is 
due to our convention for the velocity variables Vi (see 



Sec. Ill F3 1 and is explained as follows. When the particle 
collides with the piston, vi n is given in the lab frame, 
whereas v out refers to the velocity in the piston frame, 
therefore v ou t = Vi n — v p . For the subsequent bounce 
against the fixed wall, Vi n is given in the piston frame, but 
v ou t refers to the velocity in the lab frame, hence again 

Vout Vi 



Vp if the collision is elastic. Overall the distri- 



bution Eq. ( 33 1 is valid for bounces on the wall and on the 
piston side. However, when the initial velocity is nega- 
tive, the first collision takes place with the fixed wall such 
that both incoming and out-going velocities refer to the 
lab frame. In this case only, pb(v2) = 8{v2 — v\). Using 
these distributions, the bouncing times tk({vi, x) 
in Eq. ( 19 1 are reduced to the correct elastic counterpart 



tk(vi,x)and the averages, Eqs. (23 1 and (p5|, are evalu 



ated by integrating over the initial velocity v\ and initial 
position x. In this case the non-equilibrium work the- 
orem has been proven to hold exactly for all parameter 
values [BJ. 

We observe a particularly interesting feature of the 
non-equilibrium work theorem for the adiabatic expan- 
sion of the ideal gas, that is worth highlighting. In the 
quasistatic limit of a very slow moving piston, the work 
average (W) yields the free energy difference of the adi- 
abatic expansion of the gas as one would expect. The 
exponential work average, on the other hand, yields the 
free energy difference of the isothermal expansion of the 
gas. This can clearly be seen in Fig. [7] 

The averages — log (exp(— W)) and (W) respectively 
show excellent agreement with the isothermal free energy 
difference Eq. ( 18 ) (solid line) and the adiabatic free en- 
ergy difference (dotted line). The latter reads explicitly: 



1 



L 



ad 



L + V p T 



(34) 



Thus, by performing an adiabatic experiment, we ob- 
tain information about both an adiabatic and an isother- 
mal system. If we consider the irreversible case for this 
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FIG. 7: Reversible adiabatic expansion of the ideal gas for 
v p = 0.01. Both the free energy difference for isothermal 
expansion and that for adiabatic expansion are obtained from 
the same work values. 



— AF isothermal 
— - AF adiabatic 
v -ln<exp(-w)> 

o <w> 




FIG. 8: Irreversible adiabatic expansion of the ideal gas for 
v p = 1. The dissipated work is observed. 



V. CONCLUSION 

The main result of our investigation is the validation, 
both numerically and analytically, of the non-equilibrium 
work theorem for the isothermal expansion of an ideal gas 
against a piston. Although the analytical calculation is 
restricted to the limit of a fast-moving piston and small 
volume extension, the simulation confirms the result for 
a wide range of parameter values. The two main char- 
acteristics of the model under consideration should be 
emphasized again. First, the isothermal model exhibits 
strong thermal coupling between system and heat reser- 
voir, an important and more physically relevant exten- 
sion to the ideal gas models previously discussed in the 
literature [HIE]. Second, both the system and the heat 
reservoir depend on the work parameter A, violating the 
assumptions in Jarzynski's original derivation 8 . Fur- 
thermore it has been shown that microscopic reversibility 
is broken due to the moving and thermostatting piston, 
such that Crooks' derivation [9] (which assumes transi- 
tion rates independent of A) does not hold either. We 
have thus identified a regime where Jarzynski's relation 
might have been expected to fail, however it appears that 
the validity of the non-equilibrium work theorem is not 
affected. 

It is a pleasure to acknowledge helpful discussions with 
C. Jarzynski and comments by A. Adib. This work 
was funded by EPSRC Grant GR/T24593/01. RMLE 
is funded by the Royal Society. 



APPENDIX A: THE n = 1 APPROXIMATION 

We show the convergence of the averages (e°) , + 
(e°)_ — > 1 and (e~ W2 )_ — > in the limit of a fast moving 
piston and small extended volume. In order to calculate 
the averages, we resolve the theta functions with respect 
to V\ and obtain multiple dependent integrals. 

First the zero-work contribution for positive initial ve- 
locity: 



1 



+ 



2tt£ J 



L p(L~x)/T+v p 

dx / dvi e^/ 2 . (Al) 

Jo 



adiabatic process, we observe again the dissipated work 
Wd (see Fig. [8]) which is responsible for the deviation of 
(W) from AF a j. On the other hand AF (the isother- 
mal result) can still be determined by evaluating the 
exponential work average. By application of the non- 
equilibrium work theorem we can perform an adiabatic 
or an isothermal experiment and obtain the same result, 
the isothermal free energy difference. From a numerical 
point of view the isothermal simulation proves slightly 
more advantageous because more sampling takes place 
during each realization of the protocol leading to a faster 
convergence of the exponential average. 



The zero-work contribution for negative initial velocity 
reads 



1 



L />oo 

dx dv 2 v 2 e~ v */ 2 

2nL Jo J l/t+v p 

r(v 2 x/(v2-v p )T-L) 

la 

1 j-L/T+Vp 

dv 2 v 2 e""^ 2 , (A2) 





and the one-bounce contribution for negative initial ve- 
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locity is given as 



<e-«*>_ 



dx 



^/2ttL Jq Jl/t+v 

V2x/((v 2 —v p )t—L) 



dv 2 V 2 e -(^-"p) 2 /2 



x / dv 3 v 3 e 
Jo 



-(i-3-"p) 2 /2 



(A3) 



We discuss the limiting behaviour as follows. The in- 
tegration over v\ always leads to the error function 
(Eq. (Al I and Eq. (A2|) or to a sum of a constant and 
the error function (Eq. (A3)). These expressions are al- 



ways bounded by a constant for arbitrary values of the 
For Eq. (Al) we immediately obtain, in the 

< e °> + 



argument 

limit of small t and 



arge v p : {e") + -> 1/2. 

For the two other cases we see that the integrand of 
the second integration is multiplied by a bounded term. 
Thus the convergence of both averages depends only on 
the ^-integration. If we note that xexp(— x 2 ) is already 



small for x > 1, we see that the first term in Eq. ( A2 1 is 



clearly vanishing in the considere d lim it and the result is 
(e°)_ — * 1/2. In the case of Eq. (A3 1 we simply observe 
that there is an additional shift v 2 — v p in the argument 
of the exponential such that this integrand decays to zero 
even faster than the first term in Eq. ( |A2| . As a result 
the relations Eq. ( [30| are valid and the single bounce 
approximation analytically confirms the non-equilibrium 
work-theorem in the regime L v p t 3> t. 



Q(v a - Ax/At)Q(v b - Ax/ At) 
xQ(Ax/v a + Ax/v b - At)[Ax - (At - Ax/v b )v a ] 
+S(v a - Ax / Ai)Q(Ax / At - v b )Ax 
+B(Ax/At - v a )0(v b - Ax/At)Axv a /v b 
+Q(Ax/At - v a )e(Ax/At - v b )Atv a 
B(v a ,v b ;Ax,At). (Bl) 



The joint probability p(b, At; a, 0) then reads: 



p(b,At;a,0) = p eq (x, v a )p B (v b )Av 2 B(v a , v b ; Ax, At). 



(B2) 



Dividing by p(a, 0) = p eq (x, v a ) AxAv we thus obtain the 
conditional probability p(b, At\a, 0): 



p(b,At\a,0) = p B (v b )^B(v a ,v b ;Ax,At). (B3) 
Ax 



Finally this can be rewritten in the form of Eq. (13 1 if 



we consider that psivb) — v b exp[— v 2 /2] (for k B T set to 
unity) and 



APPENDIX B: THE CONDITIONAL 
PROBABILITY p(b, At\a, 0) 

The integration over the product of the two theta func- 
tions in Eq. (12 1 yields the result: 



dx 9(At - x/v a )Q(Ax - (At - x/v a )v b ) 



Av 

A(v a ,v b ) = —v b B(v a ,v b ;Ax,At). 
Ax 



(B4) 



From Eq. (Bl I we see that v b B(v a , v b ; Ax, At) is invariant 



under the exchange v a <-> v b . 
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